POGIL data analysis (all multilevel models with race and gender)
Sean Raleigh
April 17, 2023
Preliminaries
Load packages.
library(tidyverse)
library(rstan)
library(tidybayes)
library(patchwork)Set Stan options:
options(mc.cores = parallel::detectCores()) # parallel proc
rstan_options(auto_write = TRUE) # avoid recompilationImport data
data_raw <- read_csv("data/2023-01-21.csv")Rows: 1137 Columns: 10── Column specification ──────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (7): term, phase, filename, course, race, gender, grade
dbl (3): cohort, faculty, count
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Wrangle data
Remove unneeded variables.
d <- data_raw %>%
select(-filename, -term, -cohort, -course) %>%
arrange(faculty)
d$faculty <- as_factor(d$faculty)
dd %>%
group_by(faculty) %>%
count()Duplicate rows
Duplicate rows according to count variable to create one
row per student.
d_all <- d %>%
uncount(count)
d_allCheck that there are the correct number of rows:
sum(d$count)[1] 3632NROW(d_all)[1] 3632Drop missing race data
Some of the race data is missing, so we have to drop those rows.
NROW(d_all)[1] 3632d_all <- d_all %>%
drop_na()
NROW(d_all)[1] 3416Combine phases
d_all$phase_coll <- d_all$phase %>%
fct_collapse("1" = c("1a", "1b"), "2" = "2a", "3" = c("3a", "3b"))
table(d_all$phase)
1a 1b 2a 3a 3b
660 863 1024 529 340 table(d_all$phase_coll)
1 2 3
1523 1024 869 Collapse “failing” grades
d_all$grade_coll <- d_all$grade %>%
fct_collapse(DFW = c("D", "F", "W"))
table(d_all$grade)
A B C D F W
1095 829 538 254 298 402 table(d_all$grade_coll)
A B C DFW
1095 829 538 954 We’ll also need this variable as a sequence of numbers:
d_all <- d_all %>%
mutate(grade_ord = as.integer(grade_coll))
table(d_all$grade_ord)
1 2 3 4
1095 829 538 954 Collapse race variable
d_all <- d_all %>%
mutate(race_coll = ifelse(race == "nonURM", "White or Asian", "Black, Hispanic, or Indigenous"))
table(d_all$race)
amerind black hispanic natpac nonURM
17 514 429 13 2443 table(d_all$race_coll)
Black, Hispanic, or Indigenous White or Asian
973 2443 Create index variables for phases, faculty, gender, and race
Phase
The phase indicator has to be 1 or 2. Phase 1 will always be 1, but in the two different models, the comparison phase is either Phase 2 or Phase 3. Both will need to be labeled with the index 2.
d_all <- d_all %>%
mutate(phase_ind = ifelse(phase_coll == 1, 1, 2))
table(d_all$phase_coll)
1 2 3
1523 1024 869 table(d_all$phase_ind)
1 2
1523 1893 Faculty
Currently, the faculty are numbered as follows:
table(d_all$faculty)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
81 156 122 171 92 224 326 155 180 690 100 47 48 34 84 80 19 37 41 36
21 22 23 24 25
39 89 252 264 49 The faculty is a factor variable, but we’ll also need an
integer version. We’ll also drop the unused levels (faculty with no race
data at all).
table(d_all$faculty)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
81 156 122 171 92 224 326 155 180 690 100 47 48 34 84 80 19 37 41 36
21 22 23 24 25
39 89 252 264 49 d_all$faculty <- fct_drop(d_all$faculty)
d_all <- d_all %>%
mutate(faculty_ind = as.integer(as_factor(faculty)))
table(d_all$faculty_ind)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
81 156 122 171 92 224 326 155 180 690 100 47 48 34 84 80 19 37 41 36
21 22 23 24 25
39 89 252 264 49 Gender
levels(as_factor(d_all$gender))[1] "men" "women"So in the index variable, 1 = “men” and 2 = “women”.
d_all <- d_all %>%
mutate(gender_ind = as.integer(as_factor(gender)))
table(d_all$gender)
men women
2561 855 table(d_all$gender_ind)
1 2
2561 855 Race
levels(as_factor(d_all$race_coll))[1] "White or Asian" "Black, Hispanic, or Indigenous"So in the index variable, 1 = “White or Asian”, 2 = “Black, Hispanic, or Indigenous”.
d_all <- d_all %>%
mutate(race_ind = as.integer(as_factor(race_coll)))
table(d_all$race_coll)
Black, Hispanic, or Indigenous White or Asian
973 2443 table(d_all$race_ind)
1 2
2443 973 Interaction terms
Phase_gender
The new phase_gender variable is indexed as follows:
- Phase 1, Men
- Phase 1, Women
- Phase 2/3, Men
- Phase 2/3, Women
d_all <- d_all %>%
mutate(phase_gender = case_when(phase_ind == 1 & gender_ind == 1 ~ 1,
phase_ind == 1 & gender_ind == 2 ~ 2,
phase_ind == 2 & gender_ind == 1 ~ 3,
phase_ind == 2 & gender_ind == 2 ~ 4))
table(d_all$phase_ind, d_all$gender)
men women
1 1164 359
2 1397 496table(d_all$phase_gender)
1 2 3 4
1164 359 1397 496 Phase_race
The new phase_race variable is indexed as follows:
- Phase 1, White or Asian
- Phase 1, Black, Hispanic, or Indigenous
- Phase 2/3, White or Asian
- Phase 2/3, Black, Hispanic, or Indigenous
d_all <- d_all %>%
mutate(phase_race = case_when(phase_ind == 1 & race_ind == 1 ~ 1,
phase_ind == 1 & race_ind == 2 ~ 2,
phase_ind == 2 & race_ind == 1 ~ 3,
phase_ind == 2 & race_ind == 2 ~ 4))
table(d_all$phase_ind, d_all$race_coll)
Black, Hispanic, or Indigenous White or Asian
1 398 1125
2 575 1318table(d_all$phase_race)
1 2 3 4
1125 398 1318 575 Complete data set
d_allExploratory data analysis
Gender
ggplot(d_all, aes(x = gender)) +
geom_bar()
Race
ggplot(d_all, aes(x = race)) +
geom_bar()
ggplot(d_all, aes(x = race_coll)) +
geom_bar()
Grade distributions
All letter grades
All grades given across the entire study:
ggplot(d_all, aes(x = grade)) +
geom_bar()
Overall grade distributions for each phase
ggplot(d_all, aes(x = grade)) +
geom_bar() +
facet_grid(phase ~ .)
As proportions within each phase:
d_all %>%
group_by(phase, grade) %>%
count() %>%
group_by(phase) %>%
mutate(prop = n/sum(n)) %>%
ggplot(aes(x = grade, y = prop)) +
geom_bar(stat = "identity") +
facet_grid(phase ~ .)
Grade distributions by faculty across phases
ggplot(d_all, aes(x = grade)) +
geom_bar() +
facet_grid(phase ~ faculty)
As proportions:
d_all %>%
group_by(phase, faculty, grade) %>%
count() %>%
group_by(phase, faculty) %>%
mutate(prop = n/sum(n)) %>%
ggplot(aes(x = grade, y = prop)) +
geom_bar(stat = "identity") +
facet_grid(phase ~ faculty)
Grade distributions by gender
ggplot(d_all, aes(x = grade_coll)) +
geom_bar() +
facet_grid(phase_coll ~ gender)
As proportions:
d_all %>%
group_by(phase_coll, gender, grade_coll) %>%
count() %>%
group_by(phase_coll, gender) %>%
mutate(prop = n/sum(n)) %>%
ggplot(aes(x = grade_coll, y = prop)) +
geom_bar(stat = "identity") +
facet_grid(phase_coll ~ gender)
Check the numbers:
Phase 1 to Phase 2:
d_all %>%
filter(phase_coll == 1 | phase_coll == 2) %>%
group_by(grade_coll, gender, phase_coll) %>%
count() %>%
group_by(phase_coll, gender) %>%
mutate(prop = n/sum(n))d_all %>%
filter(phase_coll == 1 | phase_coll == 2) %>%
group_by(grade_coll, gender, phase_coll) %>%
count() %>%
group_by(phase_coll, gender) %>%
mutate(prop = n/sum(n)) %>%
ggplot(aes(y = prop, x = phase_coll)) +
geom_point() +
facet_grid(gender ~ grade_coll)
Phase 1 to Phase 3:
d_all %>%
filter(phase_coll == 1 | phase_coll == 3) %>%
group_by(grade_coll, gender, phase_coll) %>%
count() %>%
group_by(phase_coll, gender) %>%
mutate(prop = n/sum(n))d_all %>%
filter(phase_coll == 1 | phase_coll == 3) %>%
group_by(grade_coll, gender, phase_coll) %>%
count() %>%
group_by(phase_coll, gender) %>%
mutate(prop = n/sum(n)) %>%
ggplot(aes(y = prop, x = phase_coll)) +
geom_point() +
facet_grid(gender ~ grade_coll)
Grade distributions by race
ggplot(d_all, aes(x = grade_coll)) +
geom_bar() +
facet_grid(phase_coll ~ race_coll)
As proportions:
d_all %>%
group_by(phase_coll, race_coll, grade_coll) %>%
count() %>%
group_by(phase_coll, race_coll) %>%
mutate(prop = n/sum(n)) %>%
ggplot(aes(x = grade_coll, y = prop)) +
geom_bar(stat = "identity") +
facet_grid(phase_coll ~ race_coll)
Check the numbers:
Phase 1 to Phase 2:
d_all %>%
filter(phase_coll == 1 | phase_coll == 2) %>%
group_by(grade_coll, race_coll, phase_coll) %>%
count() %>%
group_by(phase_coll, race_coll) %>%
mutate(prop = n/sum(n))d_all %>%
filter(phase_coll == 1 | phase_coll == 2) %>%
group_by(grade_coll, race_coll, phase_coll) %>%
count() %>%
group_by(phase_coll, race_coll) %>%
mutate(prop = n/sum(n)) %>%
ggplot(aes(y = prop, x = phase_coll)) +
geom_point() +
facet_grid(race_coll ~ grade_coll)
Phase 1 to Phase 3:
d_all %>%
filter(phase_coll == 1 | phase_coll == 3) %>%
group_by(grade_coll, race_coll, phase_coll) %>%
count() %>%
group_by(phase_coll, race_coll) %>%
mutate(prop = n/sum(n))d_all %>%
filter(phase_coll == 1 | phase_coll == 3) %>%
group_by(grade_coll, race_coll, phase_coll) %>%
count() %>%
group_by(phase_coll, race_coll) %>%
mutate(prop = n/sum(n)) %>%
ggplot(aes(y = prop, x = phase_coll)) +
geom_point() +
facet_grid(race_coll ~ grade_coll)
Filter data for phase analysis
# Phase 1 and Phase 2 data
d_final_12 <- d_all %>%
filter(as.integer(phase_coll) == 1 | as.integer(phase_coll) == 2)
# Only faculty labeled 1-11 have Phase 3 data
d_final_13 <- d_all %>%
filter(as.integer(phase_coll) == 1 | as.integer(phase_coll) == 3) %>%
filter(as.integer(faculty) >= 1 & as.integer(faculty) <= 11)
Drop unused levels of phase_coll and
faculty:
d_final_12$phase_coll <- fct_drop(d_final_12$phase_coll)
d_final_13$phase_coll <- fct_drop(d_final_13$phase_coll)
d_final_12$faculty <- fct_drop(d_final_12$faculty)
d_final_13$faculty <- fct_drop(d_final_13$faculty)Final data sets
d_final_12str(d_final_12)tibble [2,547 × 15] (S3: tbl_df/tbl/data.frame)
$ phase : chr [1:2547] "1a" "1a" "1a" "1a" ...
$ faculty : Factor w/ 25 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
$ race : chr [1:2547] "nonURM" "nonURM" "nonURM" "nonURM" ...
$ gender : chr [1:2547] "men" "men" "men" "men" ...
$ grade : chr [1:2547] "A" "A" "A" "A" ...
$ phase_coll : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
$ grade_coll : Factor w/ 4 levels "A","B","C","DFW": 1 1 1 1 1 1 1 1 2 2 ...
$ grade_ord : int [1:2547] 1 1 1 1 1 1 1 1 2 2 ...
$ race_coll : chr [1:2547] "White or Asian" "White or Asian" "White or Asian" "White or Asian" ...
$ phase_ind : num [1:2547] 1 1 1 1 1 1 1 1 1 1 ...
$ faculty_ind : int [1:2547] 1 1 1 1 1 1 1 1 1 1 ...
$ gender_ind : int [1:2547] 1 1 1 1 1 2 2 2 1 1 ...
$ race_ind : int [1:2547] 1 1 1 1 1 1 1 1 1 1 ...
$ phase_gender: num [1:2547] 1 1 1 1 1 2 2 2 1 1 ...
$ phase_race : num [1:2547] 1 1 1 1 1 1 1 1 1 1 ...d_final_13str(d_final_13)tibble [1,748 × 15] (S3: tbl_df/tbl/data.frame)
$ phase : chr [1:1748] "1a" "1a" "1a" "1a" ...
$ faculty : Factor w/ 11 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
$ race : chr [1:1748] "nonURM" "nonURM" "nonURM" "nonURM" ...
$ gender : chr [1:1748] "men" "men" "men" "men" ...
$ grade : chr [1:1748] "A" "A" "A" "A" ...
$ phase_coll : Factor w/ 2 levels "1","3": 1 1 1 1 1 1 1 1 1 1 ...
$ grade_coll : Factor w/ 4 levels "A","B","C","DFW": 1 1 1 1 1 1 1 1 2 2 ...
$ grade_ord : int [1:1748] 1 1 1 1 1 1 1 1 2 2 ...
$ race_coll : chr [1:1748] "White or Asian" "White or Asian" "White or Asian" "White or Asian" ...
$ phase_ind : num [1:1748] 1 1 1 1 1 1 1 1 1 1 ...
$ faculty_ind : int [1:1748] 1 1 1 1 1 1 1 1 1 1 ...
$ gender_ind : int [1:1748] 1 1 1 1 1 2 2 2 1 1 ...
$ race_ind : int [1:1748] 1 1 1 1 1 1 1 1 1 1 ...
$ phase_gender: num [1:1748] 1 1 1 1 1 2 2 2 1 1 ...
$ phase_race : num [1:1748] 1 1 1 1 1 1 1 1 1 1 ...Final summary tables
table(d_all$gender, d_all$race)
amerind black hispanic natpac nonURM
men 14 366 321 11 1849
women 3 148 108 2 594table(d_final_12$race_coll, d_final_12$gender, d_final_12$phase_ind), , = 1
men women
Black, Hispanic, or Indigenous 294 104
White or Asian 870 255
, , = 2
men women
Black, Hispanic, or Indigenous 210 81
White or Asian 542 191table(d_final_12$race_coll, d_final_12$gender, d_final_12$phase_ind) %>%
prop.table(margin = 3), , = 1
men women
Black, Hispanic, or Indigenous 0.19304005 0.06828628
White or Asian 0.57124097 0.16743270
, , = 2
men women
Black, Hispanic, or Indigenous 0.20507812 0.07910156
White or Asian 0.52929688 0.18652344table(d_final_13$race_coll, d_final_13$gender, d_final_13$phase_ind), , = 1
men women
Black, Hispanic, or Indigenous 170 54
White or Asian 514 141
, , = 2
men women
Black, Hispanic, or Indigenous 208 76
White or Asian 437 148table(d_final_13$race_coll, d_final_13$gender, d_final_13$phase_ind) %>%
prop.table(margin = 3), , = 1
men women
Black, Hispanic, or Indigenous 0.19340159 0.06143345
White or Asian 0.58475540 0.16040956
, , = 2
men women
Black, Hispanic, or Indigenous 0.23935558 0.08745685
White or Asian 0.50287687 0.17031070Models
The ordinal regression model assumes a continuous latent distribution for the response variable with “cutpoints” that generate ordinal categories. For grades, this makes a lot of sense, as the continuous distribution of points is converted into discrete ordered grades.
Suppose there are \(K > 2\) categories in our ordinal response \(Y\). (For the grade data, \(K = 4\).) Let \(c\) be a vector of length \(K + 1\) with \(c_{0} = -\infty\), \(c_{K} = \infty\), and \(c_{1} < c_{2} < \dots < c_{K-1}\). Finally, let \(\eta\) be a real number. Then for \(1 \leq k \leq K\), the ordered logistic distribution is given by
\[ OrdLog(k \mid \eta, c) = invlogit(\eta - c_{k - 1}) - invlogit(\eta - c_{k}) \] where \(invlogit\) is the inverse logit function (sometimes called the logistic function):
\[ invlogit(x) = \frac{e^{x}}{1 + e^{x}}. \]
Note that \(invlogit(-\infty) = 0\) and \(invlogit(\infty) = 1\).
\(OrdLog(k \mid \eta, c)\) is the probability (given the model) of the ordinal response variable taking the value \(k\): \(Pr(Y = k)\). (Although \(Y\) consists of ordered categories, these categories are assigned integers for convenience in writing down the model mathematically. The model does not assume that the categories are equally spaced.)
This is not the most natural parameterization of the ordinal regression model, but it’s the one Stan uses, so we’re stuck with it. We can re-express the model other ways that are easier to understand. First, note that the inverse logit function represents the cumulative probability of a logistic distributed random variable \(X\); in other words,
\[ invlogit(x) = Pr(X \leq x) \] where \[ X \sim Logistic(0, 1). \]
Replacing \(x\) with \(\eta - c_{k-1}\) and using the identity \(invlogit(-x) = 1 - invlogit(x)\):
\[\begin{align*} Pr(Y = k) &= invlogit(\eta - c_{k - 1}) - invlogit(\eta - c_{k}) \\ &= \left( 1 - invlogit(c_{k - 1} - \eta) \right) - \left( 1 - invlogit(c_{k} - \eta) \right) \\ &= invlogit(c_{k} - \eta) - invlogit(c_{k - 1} - \eta) \end{align*}\]
One can show (either using algebra and the expression above, or thinking about \(invlogit\) as a cumulative probability) that
\[ Pr(Y \leq k) = invlogit(c_{k} - \eta) \]
Taking the logit function (log odds) on both sides, we get an equivalent expression representing a linear model predicting the cumulative log odds of the ordinal response, meaning the log odds of being in category \(k\) or lower.
\[ logit(Pr(Y \leq k)) = \log\left( \frac{Pr(Y \leq k)}{Pr(Y > k)} \right) = c_{k} - \eta \]
Based on the way grades are parameterized, however, lower values of \(k\) are actually better grades. For example, the event \(Y \leq 2\) refers to grades of B or A.
The effect of \(\eta\) is as follows: if \(\eta > 0\), then the log odds—and, therefore, the probability—of being in category \(k\) or lower (meaning higher grades) will be smaller. This shifts probability to higher-numbered categories. In other words, larger values of \(\eta\) will correspond to more probability assigned to higher-numbered categories of the response (corresponding to lower grades). Likewise, negative values of \(\eta\) will shift more probability onto the lower-numbered categories (higher grades). One important point to make is that \(\eta\) does not depend on \(k\). The “shift” in log-odds is assumed by the model to be constant across all values of \(k\). This is called the “proportional odds” assumption due to the following equation:
\[ \frac{Pr(Y \leq k)}{Pr(Y > k)} = e^{(c_{k} - \eta)} = e^{c_{k}} e^{\eta} \]
In general, predictors in the model will be incorporated into the \(\eta\) term. In this data, we have the following parameters coded as index variables (one variable for each group in the data):
- \(\alpha_{faculty}\): Each faculty participant has a parameter. These are so called “random effects” and are modeled hierarchically using a normal hyperprior: \(\alpha \sim Normal(\bar{\alpha}, \sigma)\). The result will be a partial pooling of the results across faculty.
- \(\beta_{phase}\): Each phase has a parameter. There are two phases. In one model, they are Phase 1 and Phase 2. In the other model, they are Phase 1 and Phase 3.
- \(\gamma_{gender}\): Genders are “men” and “women”.
- \(\rho_{race}\): Racial categories are “Black, Hispanic, or Indigenous” and “White or Asian”.
- \(\delta_{phase/gender}\): These are interaction terms that model deviations from the main effects for phase and gender. (There are four such terms, one for each combination of phase and gender.)
- \(\zeta_{phase/race}\): These are interaction terms that model deviations from the main effects for phase and race. (There are four such terms, one for each combination of phase and race.)
There are two models we are considering. The first model (Model A) uses
\[ \eta = \alpha_{faculty} + \beta_{phase} \]
This model will measure the average effect between phases and will allow us to summarize inferences about grade changes for individual faculty and for all faculty collectively.
The second model (Model B) will consider all predictors, incorporating race and gender, including interaction terms with phase:
\[ \eta = \alpha_{faculty} + \beta_{phase} + \gamma_{gender} + \rho_{race} + \delta_{phase/gender} + \zeta_{phase/race} \]
Model B will measure differences between gender and race categories. In theory, it should be possible to recover Model A from Model B by marginalizing over gender and race. However, the parameterization of the interaction terms in Model B makes it computationally difficult to perform that marginalization. Therefore, we opt to run both Model A and Model B.
Suppose that \(Y_{1}\) is the ordinal response (grades) for Phase 1 and \(Y_{j}\) is the ordinal response for either Phase 2 or Phase 3 (\(j = 2 \text{ or } 3\)). For each cutpoint \(c_{k}\):
\[ \eta = c_{k} - logit(Pr(Y \leq k)). \]
In Model A, we then have the following equations (for all faculty):
\[\begin{align*} \beta_{1} &= c_{k} - \alpha_{faculty} - logit(Pr(Y_{1} \leq k)) \\ \beta_{j} &= c_{k} - \alpha_{faculty} - logit(Pr(Y_{j} \leq k)) \end{align*}\]
We can subtract these two equations to obtain
\[\begin{align*} \beta_{1} - \beta_{j} &= logit(Pr(Y_{j} \leq k)) - logit(Pr(Y_{1} \leq k)) \\ &= log\left( \frac{Pr(Y_{j} \leq k)}{Pr(Y_{j} > k)} \right) - log\left( \frac{Pr(Y_{1} \leq k)}{Pr(Y_{1} > k)} \right) \\ &= log\left( \frac{Pr(Y_{j} \leq k)/Pr(Y_{j} > k)}{Pr(Y_{1} \leq k)/Pr(Y_{1} > k)} \right) \end{align*}\]
Therefore, we can exponentiate the expression on the left to get the odds ratio from Phase 1 to Phase 2/3.
An odds ratio of 1 indicates that the grade distribution does not change from one phase to the next. An odds ratio greater than 1 indicates that grades improved from Phase 1 to Phase \(j\). An odds ratio less than 1 indicates that grades decreased from Phase 1 to Phase \(j\). The proportional odds assumption forces this odds ratio to be the same at all grade levels. Although this is rarely a realistic assumption, the model effectively averages over any differences from one grade level to the next to enforce this assumption.
In Model B, we have the following equations for each value of faculty, gender, and race:
\[\begin{align*} \beta_{1} + \delta_{1/gender} + \zeta_{1/race} &= c_{k} - \alpha_{faculty} - \gamma_{gender} - \rho_{race} - logit(Pr(Y_{1} \leq k)) \\ \beta_{j} + \delta_{j/gender} + \zeta_{j/race} &= c_{k} - \alpha_{faculty} - \gamma_{gender} - \rho_{race} - logit(Pr(Y_{j} \leq k)) \end{align*}\]
We can subtract these two equations to obtain
\[\begin{align*} \beta_{1} - \beta_{j} + \delta_{1/gender} - \delta_{j/gender} + \zeta_{1/race} - \zeta_{j/race} &= logit(Pr(Y_{j} \leq k)) - logit(Pr(Y_{1} \leq k)) \\ &= log\left( \frac{Pr(Y_{j} \leq k)}{Pr(Y_{j} > k)} \right) - log\left( \frac{Pr(Y_{1} \leq k)}{Pr(Y_{1} > k)} \right) \\ &= log\left( \frac{Pr(Y_{j} \leq k)/Pr(Y_{j} > k)}{Pr(Y_{1} \leq k)/Pr(Y_{1} > k)} \right) \end{align*}\]
Just as before, we can exponentiate the expression on the left to get the odds ratio from Phase 1 to Phase 2/3 for but now for any combination of gender and race.
As this is a Bayesian analysis, we must specify priors for all parameters. We choose fairly generic, weakly-informative priors. As they are centered at 0, they will also serve as regularizing priors. (The data is sufficient in this study that they will overwhelm just about any sensible choice of prior.) The priors are as follows:
\[\begin{align*} \alpha_{faculty} &\sim Normal(\bar{\alpha}, \sigma) \\ \bar{\alpha} &\sim Normal(0, 1) \\ \sigma &\sim Exponential(1) \\ \beta_{phase} &\sim Normal(0, 1) \\ \gamma_{gender} &\sim Normal(0, 1) \\ \rho_{race} &\sim Normal(0, 1) \\ \delta_{phase/gender} &\sim Normal(0, 1) \\ \zeta_{phase/race} &\sim Normal(0, 1) \\ \end{align*}\]
The components of the cutpoint vector are each assigned a uniform prior by Stan (which is truncated based on the constraints in the ordering).
Multilevel models in Stan
Model A
data {
int<lower=0> N; // sample size
int<lower=1, upper=2> phase_ind[N]; // two phases
int<lower=2> J; // number of faculty
int<lower=1> faculty_ind[N]; // faculty array
int<lower=2> K; // number of grade categ.
int<lower=1, upper=K> y[N]; // grade codes
}
parameters {
real alpha_bar;
real<lower=0> sigma;
real alpha[J];
real beta[2];
ordered[K - 1] c;
}
model {
alpha_bar ~ normal(0, 1);
sigma ~ exponential(1);
alpha ~ normal(alpha_bar, sigma);
beta ~ normal(0, 1);
for (n in 1:N) {
real eta;
eta = alpha[faculty_ind[n]] + beta[phase_ind[n]];
y[n] ~ ordered_logistic(eta, c);
}
}Phase 1 to Phase 2
The data list:
ord_mm_dataA_12 <- list(
N = NROW(d_final_12),
phase_ind = d_final_12$phase_ind,
J = nlevels(d_final_12$faculty),
faculty_ind = d_final_12$faculty_ind,
K = nlevels(d_final_12$grade_coll),
y = d_final_12$grade_ord
)Fit the model:
ord_mm_fitA_12 <- sampling(ord_mm_modelA,
data = ord_mm_dataA_12,
iter = 8000,
seed = 11114)Coefficients of the model:
summary(ord_mm_fitA_12, pars = c("alpha_bar", "sigma", "alpha", "beta"))$summary mean se_mean sd 2.5% 25% 50%
alpha_bar 0.029001437 0.029786871 0.9536399 -1.8069434 -0.61646407 0.03504065
sigma 0.670315464 0.001652577 0.1135559 0.4842918 0.59108139 0.65850768
alpha[1] 0.075343637 0.030426792 0.9847952 -1.8417903 -0.58778038 0.07961142
alpha[2] -0.102168317 0.030316231 0.9756223 -1.9931740 -0.75341640 -0.10033212
alpha[3] 0.437905864 0.030466097 0.9829083 -1.4700580 -0.21648871 0.44384830
alpha[4] 1.517049244 0.030286087 0.9773438 -0.3752865 0.85733245 1.51134438
alpha[5] 0.834573605 0.030410173 0.9814089 -1.0883020 0.16892742 0.83439434
alpha[6] 0.591165648 0.030373492 0.9747921 -1.2891085 -0.06682071 0.59204394
alpha[7] 0.240242682 0.030385629 0.9695838 -1.6472710 -0.40941128 0.23948115
alpha[8] 0.476283274 0.030259333 0.9762928 -1.4037574 -0.16916571 0.47107418
alpha[9] -0.412525274 0.030381400 0.9836608 -2.3279657 -1.06832987 -0.40997261
alpha[10] -0.009407759 0.030430601 0.9687638 -1.8676814 -0.65811327 -0.01403333
alpha[11] -0.952224879 0.030456435 0.9841136 -2.8609631 -1.61428939 -0.94908875
alpha[12] 0.112525429 0.030542971 0.9908766 -1.7959393 -0.55006481 0.12324985
alpha[13] -0.484842878 0.030383686 0.9897050 -2.3869730 -1.15515869 -0.48206765
alpha[14] -0.828724981 0.030542724 1.0005741 -2.7864126 -1.50233594 -0.81760096
alpha[15] -0.872364325 0.030515755 0.9811754 -2.7695503 -1.53546168 -0.86629816
alpha[16] 0.395148834 0.030530467 0.9840926 -1.5133611 -0.26054166 0.39675461
alpha[17] -0.270767884 0.030604498 1.0259388 -2.2828263 -0.95767374 -0.27159709
alpha[18] -0.699694944 0.030559616 0.9975773 -2.6256240 -1.36684413 -0.69544010
alpha[19] 0.275827467 0.030482716 0.9898091 -1.6488446 -0.39167382 0.28635978
alpha[20] -1.016296672 0.030472548 1.0061288 -2.9940627 -1.67855882 -1.01097200
alpha[21] 0.227325127 0.030447877 0.9963309 -1.7179675 -0.44275006 0.23506741
alpha[22] 0.051852208 0.030416069 0.9794981 -1.8508505 -0.61139965 0.05423709
alpha[23] 0.406150652 0.030284014 0.9677748 -1.4687780 -0.24031430 0.40669707
alpha[24] 0.351915456 0.030408852 0.9702792 -1.5216778 -0.30141692 0.35464972
alpha[25] 0.376041968 0.030571050 0.9880088 -1.5098196 -0.28862058 0.37792905
beta[1] -0.020511886 0.012419711 0.7177639 -1.4185088 -0.50887378 -0.02490622
beta[2] -0.012251241 0.012431407 0.7180507 -1.4054210 -0.49855818 -0.01805745
75% 97.5% n_eff Rhat
alpha_bar 0.65932099 1.9359699 1024.989 1.003871
sigma 0.73662522 0.9264071 4721.677 1.000199
alpha[1] 0.73274718 2.0634933 1047.561 1.004016
alpha[2] 0.54846165 1.8396640 1035.650 1.003887
alpha[3] 1.08963351 2.3857676 1040.860 1.004041
alpha[4] 2.17290742 3.4643214 1041.378 1.003650
alpha[5] 1.48332315 2.7759597 1041.507 1.003935
alpha[6] 1.23547807 2.5469370 1029.994 1.004037
alpha[7] 0.88266915 2.1991909 1018.203 1.004002
alpha[8] 1.12752455 2.4350621 1040.978 1.003686
alpha[9] 0.24415900 1.5379682 1048.275 1.003798
alpha[10] 0.63255957 1.9494239 1013.479 1.004110
alpha[11] -0.29342285 1.0265801 1044.077 1.003779
alpha[12] 0.77078097 2.1066684 1052.487 1.003858
alpha[13] 0.17139853 1.4827349 1061.037 1.003906
alpha[14] -0.16765236 1.1736633 1073.206 1.004296
alpha[15] -0.22090492 1.0806423 1033.820 1.004058
alpha[16] 1.05009349 2.3719805 1038.975 1.004060
alpha[17] 0.42311115 1.7441187 1123.757 1.004024
alpha[18] -0.04191571 1.2920292 1065.608 1.003902
alpha[19] 0.93451146 2.2354252 1054.376 1.003851
alpha[20] -0.34323283 0.9945033 1090.158 1.003453
alpha[21] 0.88491649 2.2171462 1070.763 1.003673
alpha[22] 0.70629204 2.0146864 1037.053 1.003813
alpha[23] 1.05346697 2.3513569 1021.226 1.004063
alpha[24] 0.99506489 2.3096017 1018.107 1.003973
alpha[25] 1.03320761 2.3677921 1044.482 1.004092
beta[1] 0.45326510 1.4041386 3339.952 1.001816
beta[2] 0.46502102 1.4183657 3336.335 1.001856Phase 1 to Phase 3
The data list:
ord_mm_dataA_13 <- list(
N = NROW(d_final_13),
phase_ind = d_final_13$phase_ind,
J = nlevels(d_final_13$faculty),
faculty_ind = d_final_13$faculty_ind,
K = nlevels(d_final_13$grade_coll),
y = d_final_13$grade_ord
)Fit the model:
ord_mm_fitA_13 <- sampling(ord_mm_modelA,
data = ord_mm_dataA_13,
seed = 11115)Coefficients of the model:
summary(ord_mm_fitA_13, pars = c("alpha_bar", "sigma", "alpha", "beta"))$summary mean se_mean sd 2.5% 25% 50%
alpha_bar 0.03280878 0.038745313 0.9444168 -1.8333953 -0.5908930 0.05805456
sigma 0.70510253 0.004308172 0.1853772 0.4381727 0.5727942 0.67509822
alpha[1] 0.02227800 0.040921451 0.9885992 -1.9612982 -0.6128863 0.04842078
alpha[2] -0.19695505 0.040900475 0.9799683 -2.1401574 -0.8253822 -0.18223778
alpha[3] -0.06792481 0.040916144 0.9854422 -2.0272433 -0.7191375 -0.04196462
alpha[4] 1.30440591 0.040997573 0.9765458 -0.6471014 0.6645066 1.32022014
alpha[5] 0.48140922 0.040984251 0.9882347 -1.5057686 -0.1754425 0.50338012
alpha[6] 0.43553119 0.040558890 0.9739406 -1.4894629 -0.2060175 0.45536832
alpha[7] 0.17244269 0.041279739 0.9752734 -1.7446940 -0.4567361 0.19040401
alpha[8] 0.24565720 0.040806649 0.9787358 -1.6951679 -0.3918444 0.25529683
alpha[9] -0.59432499 0.040911290 0.9812033 -2.5381528 -1.2375589 -0.58487727
alpha[10] -0.30959871 0.040891329 0.9708060 -2.2334812 -0.9395077 -0.29110418
alpha[11] -1.07363253 0.041020881 0.9927481 -3.0802238 -1.7271077 -1.07116182
beta[1] 0.11599572 0.017689802 0.7125427 -1.2554116 -0.3617500 0.12389328
beta[2] -0.11025070 0.017748871 0.7118780 -1.4700615 -0.5920231 -0.10628502
75% 97.5% n_eff Rhat
alpha_bar 0.67580689 1.8091672 594.1403 1.008116
sigma 0.80171599 1.1559964 1851.5125 1.000046
alpha[1] 0.66602255 1.9177266 583.6312 1.008871
alpha[2] 0.45407783 1.6714204 574.0734 1.008698
alpha[3] 0.59189776 1.7853584 580.0600 1.008253
alpha[4] 1.94775144 3.1623100 567.3733 1.008687
alpha[5] 1.16509885 2.3732881 581.4149 1.008136
alpha[6] 1.08499728 2.3014200 576.6241 1.008402
alpha[7] 0.82879232 2.0208206 558.1859 1.009376
alpha[8] 0.89867179 2.0997117 575.2665 1.008457
alpha[9] 0.06491474 1.2503523 575.2169 1.008848
alpha[10] 0.32764178 1.5160903 563.6409 1.008860
alpha[11] -0.41679210 0.8691978 585.6905 1.008319
beta[1] 0.59768674 1.4806983 1622.4669 1.000055
beta[2] 0.37019590 1.2717695 1608.6800 1.000321Model B
data {
int<lower=0> N; // sample size
int<lower=1, upper=2> phase_ind[N]; // two phases
int<lower=2> J; // number of faculty
int<lower=1> faculty_ind[N]; // faculty array
int<lower=1> gender_ind[N]; // gender array
int<lower=1> phase_gender[N]; // phase/gender interaction
int<lower=1> race_ind[N]; // race array
int<lower=1> phase_race[N]; // phase/race interaction
int<lower=2> K; // number of grade categories
int<lower=1, upper=K> y[N]; // grade codes
}
parameters {
real alpha_bar;
real<lower=0> sigma;
real alpha[J];
real beta[2];
real gamma[2];
real delta[4];
real rho[2];
real zeta[4];
ordered[K - 1] c;
}
model {
alpha_bar ~ normal(0, 1);
sigma ~ exponential(1);
alpha ~ normal(alpha_bar, sigma);
beta ~ normal(0, 1);
gamma ~ normal(0, 1);
delta ~ normal(0, 1);
rho ~ normal(0, 1);
zeta ~ normal(0, 1);
for (n in 1:N) {
real eta;
eta = alpha[faculty_ind[n]] + beta[phase_ind[n]] +
gamma[gender_ind[n]] + delta[phase_gender[n]] +
rho[race_ind[n]] + zeta[phase_race[n]];
y[n] ~ ordered_logistic(eta, c);
}
}Phase 1 to Phase 2
The data list:
ord_mm_dataB_12 <- list(
N = NROW(d_final_12),
phase_ind = d_final_12$phase_ind,
J = nlevels(d_final_12$faculty),
faculty_ind = d_final_12$faculty_ind,
gender_ind = d_final_12$gender_ind,
phase_gender = d_final_12$phase_gender,
race_ind = d_final_12$race_ind,
phase_race = d_final_12$phase_race,
K = nlevels(d_final_12$grade_coll),
y = d_final_12$grade_ord
)Fit the model:
ord_mm_fitB_12 <- sampling(ord_mm_modelB,
data = ord_mm_dataB_12,
iter = 8000,
seed = 11116)Coefficients of the model:
summary(ord_mm_fitB_12, pars = c("alpha_bar", "sigma", "alpha", "beta",
"gamma", "delta", "rho", "zeta"))$summary mean se_mean sd 2.5% 25%
alpha_bar -0.0268762639 0.0267624113 0.9937118 -1.9771950 -0.70317724
sigma 0.6124813049 0.0007441803 0.1076658 0.4354759 0.53674659
alpha[1] 0.1149859515 0.0272224718 1.0243387 -1.8917763 -0.58763298
alpha[2] -0.0430983029 0.0272356946 1.0130657 -2.0362678 -0.73381516
alpha[3] 0.3865731896 0.0270763595 1.0170150 -1.5986965 -0.31187190
alpha[4] 1.2810060015 0.0270861747 1.0170612 -0.7088434 0.58964716
alpha[5] 0.6496153452 0.0272142919 1.0188881 -1.3229947 -0.04826851
alpha[6] 0.3474570937 0.0271517665 1.0130124 -1.6307638 -0.34086975
alpha[7] 0.2930616644 0.0271446213 1.0087991 -1.6781631 -0.39659992
alpha[8] 0.4745274718 0.0271080374 1.0151001 -1.5022727 -0.21996119
alpha[9] -0.5416671382 0.0272853191 1.0215090 -2.5281665 -1.24147568
alpha[10] -0.1965815020 0.0271200079 1.0069573 -2.1556381 -0.88381013
alpha[11] -0.9044956803 0.0271632670 1.0216099 -2.9101445 -1.60283897
alpha[12] 0.1561437312 0.0270311130 1.0251651 -1.8345175 -0.53610502
alpha[13] -0.4682412082 0.0272591286 1.0275613 -2.4734557 -1.17032595
alpha[14] -0.7825654976 0.0272835280 1.0361773 -2.8024535 -1.48694745
alpha[15] -0.7923717001 0.0272007025 1.0195850 -2.7980248 -1.48514200
alpha[16] 0.4642855565 0.0273053707 1.0201964 -1.5282186 -0.23048393
alpha[17] -0.3001316136 0.0271014558 1.0568241 -2.3874592 -1.02456385
alpha[18] -0.7238948867 0.0273879051 1.0396533 -2.7488775 -1.43928554
alpha[19] -0.1140586466 0.0271430245 1.0267024 -2.1137006 -0.81457357
alpha[20] -0.9583003430 0.0272978352 1.0415367 -2.9937823 -1.66487007
alpha[21] 0.1982960653 0.0270619798 1.0322711 -1.8320052 -0.51086709
alpha[22] 0.0006873965 0.0272219327 1.0169526 -1.9796695 -0.70132329
alpha[23] 0.2852584204 0.0270239835 1.0077783 -1.6847825 -0.40214108
alpha[24] 0.1733491225 0.0271804376 1.0086593 -1.8073602 -0.51074801
alpha[25] 0.3218258720 0.0272212769 1.0281479 -1.6920504 -0.38839958
beta[1] -0.0365592338 0.0060725363 0.8576158 -1.7155849 -0.61248750
beta[2] 0.0397203108 0.0061355338 0.8750207 -1.7010115 -0.55906586
gamma[1] 0.0726196317 0.0063304629 0.8237138 -1.5754376 -0.47108270
gamma[2] -0.0760252516 0.0062617154 0.8188428 -1.6860356 -0.62551660
delta[1] 0.1062041681 0.0059075486 0.7863891 -1.4245766 -0.42443387
delta[2] -0.1393220636 0.0058636269 0.7845242 -1.6905555 -0.66479199
delta[3] -0.0416191068 0.0058377609 0.7743450 -1.5525482 -0.56300291
delta[4] 0.0821118475 0.0059716076 0.7844037 -1.4502454 -0.44311003
rho[1] -0.2473758000 0.0059974900 0.8015921 -1.7952746 -0.79349272
rho[2] 0.2450904501 0.0061622100 0.8107048 -1.3458526 -0.30518588
zeta[1] -0.0700634967 0.0058551843 0.7863971 -1.6128656 -0.60222550
zeta[2] 0.0298895111 0.0059422881 0.7789123 -1.5092066 -0.48981695
zeta[3] -0.1673098435 0.0058286237 0.7852595 -1.7108262 -0.69816629
zeta[4] 0.2160799118 0.0059430866 0.7788470 -1.3208207 -0.31684202
50% 75% 97.5% n_eff Rhat
alpha_bar -0.022275871 0.64182497 1.9373621 1378.702 1.0030502
sigma 0.601420957 0.67535893 0.8546775 20931.448 0.9999192
alpha[1] 0.128667966 0.80421510 2.1372176 1415.898 1.0030162
alpha[2] -0.037733402 0.63501327 1.9716809 1383.561 1.0031744
alpha[3] 0.389194738 1.06544724 2.3853421 1410.828 1.0028875
alpha[4] 1.290077798 1.96929311 3.2983117 1409.934 1.0029396
alpha[5] 0.662394959 1.32840328 2.6620387 1401.712 1.0030125
alpha[6] 0.360422325 1.02539846 2.3468869 1391.981 1.0031597
alpha[7] 0.304381573 0.97727300 2.2907300 1381.153 1.0030114
alpha[8] 0.484981447 1.16111095 2.4744262 1402.238 1.0030597
alpha[9] -0.534457622 0.14835862 1.4693843 1401.607 1.0031063
alpha[10] -0.184294946 0.48244586 1.7864202 1378.613 1.0031900
alpha[11] -0.894273613 -0.21407103 1.1239928 1414.510 1.0031428
alpha[12] 0.160934893 0.84888962 2.2006323 1438.334 1.0027778
alpha[13] -0.456218668 0.22020983 1.5531123 1420.992 1.0032547
alpha[14] -0.776918287 -0.08836516 1.2658990 1442.338 1.0028621
alpha[15] -0.787926849 -0.10117450 1.2163490 1405.033 1.0031646
alpha[16] 0.472526556 1.14939188 2.4860101 1395.955 1.0030304
alpha[17] -0.289630292 0.42376300 1.7614096 1520.618 1.0027072
alpha[18] -0.713496189 -0.03481333 1.3450132 1440.985 1.0030018
alpha[19] -0.102109514 0.57681151 1.9236175 1430.779 1.0032206
alpha[20] -0.942026281 -0.26037218 1.0783301 1455.770 1.0031342
alpha[21] 0.210445269 0.90172362 2.2389241 1455.018 1.0026026
alpha[22] 0.008686386 0.68936524 2.0096748 1395.608 1.0031388
alpha[23] 0.295872805 0.96821847 2.2761050 1390.693 1.0030078
alpha[24] 0.181276445 0.86085519 2.1628625 1377.134 1.0032437
alpha[25] 0.329772652 1.02473482 2.3584777 1426.574 1.0028032
beta[1] -0.045781289 0.53801880 1.6504596 19945.515 0.9999034
beta[2] 0.042423035 0.62880942 1.7749867 20339.111 1.0000456
gamma[1] 0.077475321 0.61973743 1.6797851 16930.966 1.0005421
gamma[2] -0.067114570 0.48076196 1.4978922 17100.721 1.0000316
delta[1] 0.105073571 0.63309246 1.6589932 17719.865 0.9999375
delta[2] -0.137477200 0.39220395 1.3933264 17901.115 0.9998856
delta[3] -0.040043546 0.47857012 1.5001161 17594.476 0.9999679
delta[4] 0.085183496 0.60695230 1.6160335 17254.274 0.9998827
rho[1] -0.254726937 0.29445746 1.3201161 17863.553 1.0002577
rho[2] 0.240647986 0.79689189 1.8474214 17308.222 0.9998542
zeta[1] -0.070088079 0.45389795 1.4614544 18038.594 0.9998077
zeta[2] 0.029407439 0.55694297 1.5524679 17181.844 1.0001003
zeta[3] -0.167745411 0.36415608 1.3621117 18150.739 0.9998427
zeta[4] 0.220400019 0.75235478 1.7127167 17174.344 1.0004454Phase 1 to Phase 3
The data list:
ord_mm_dataB_13 <- list(
N = NROW(d_final_13),
phase_ind = d_final_13$phase_ind,
J = nlevels(d_final_13$faculty),
faculty_ind = d_final_13$faculty_ind,
gender_ind = d_final_13$gender_ind,
phase_gender = d_final_13$phase_gender,
race_ind = d_final_13$race_ind,
phase_race = d_final_13$phase_race,
K = nlevels(d_final_13$grade_coll),
y = d_final_13$grade_ord
)Fit the model:
ord_mm_fitB_13 <- sampling(ord_mm_modelB,
data = ord_mm_dataB_13,
seed = 11117)Coefficients of the model:
summary(ord_mm_fitB_13, pars = c("alpha_bar", "sigma", "alpha", "beta",
"gamma", "delta", "rho", "zeta"))$summary mean se_mean sd 2.5% 25% 50%
alpha_bar -0.001471763 0.036209140 1.0090583 -2.0492627 -0.6556299 0.038849541
sigma 0.646337317 0.002902943 0.1807787 0.3854657 0.5209696 0.615585703
alpha[1] 0.167195026 0.038196086 1.0495344 -1.9227590 -0.5288647 0.190847968
alpha[2] -0.112513136 0.038418813 1.0465414 -2.2046255 -0.7960780 -0.093334432
alpha[3] -0.033576357 0.038379100 1.0492057 -2.1174890 -0.7204179 -0.007690968
alpha[4] 1.086458369 0.037872604 1.0367399 -0.9983651 0.4119324 1.107074094
alpha[5] 0.322504441 0.037954126 1.0461087 -1.8152115 -0.3627410 0.342943497
alpha[6] 0.161147217 0.037945068 1.0350511 -1.9002385 -0.5087899 0.197006419
alpha[7] 0.254616281 0.038045096 1.0337275 -1.7718010 -0.4189860 0.277420783
alpha[8] 0.284155533 0.038355120 1.0404414 -1.7911308 -0.3967445 0.316665686
alpha[9] -0.623452941 0.037530908 1.0424816 -2.7077485 -1.2889570 -0.599303291
alpha[10] -0.491427351 0.037819535 1.0306757 -2.5667558 -1.1576600 -0.475650489
alpha[11] -0.999196111 0.038256326 1.0491466 -3.0954857 -1.7093561 -0.973194244
beta[1] 0.030106392 0.011643969 0.8661377 -1.7461084 -0.5464601 0.042650300
beta[2] 0.007977114 0.011763037 0.8909157 -1.7766736 -0.5782747 0.006536847
gamma[1] 0.076963579 0.012320637 0.8423541 -1.5426931 -0.4996361 0.082187602
gamma[2] -0.088447303 0.012237844 0.8127498 -1.6847978 -0.6524259 -0.081964739
delta[1] 0.134347411 0.010490586 0.7685303 -1.3416538 -0.3911382 0.125998767
delta[2] -0.135605897 0.010199780 0.8011535 -1.6924313 -0.6871572 -0.121869196
delta[3] -0.102299909 0.010753558 0.7907185 -1.6682730 -0.6352055 -0.097133478
delta[4] 0.074091739 0.010206675 0.7822406 -1.4668425 -0.4529287 0.069059161
rho[1] -0.290383952 0.010532786 0.8113602 -1.8932760 -0.8433822 -0.285201640
rho[2] 0.305576529 0.010220241 0.8031747 -1.2708818 -0.2251190 0.305091262
zeta[1] -0.079082692 0.010237516 0.7769264 -1.5944271 -0.6095494 -0.076078534
zeta[2] 0.075866887 0.010345664 0.7688075 -1.4023087 -0.4510489 0.073055317
zeta[3] -0.229343801 0.010894157 0.7640890 -1.7390401 -0.7386338 -0.216351345
zeta[4] 0.195547966 0.010708328 0.7774079 -1.3546937 -0.3303042 0.208079063
75% 97.5% n_eff Rhat
alpha_bar 0.66547283 1.921028 776.5977 1.0038538
sigma 0.73922979 1.084112 3878.0868 0.9998272
alpha[1] 0.84713448 2.211668 755.0151 1.0040731
alpha[2] 0.58994280 1.892645 742.0360 1.0037583
alpha[3] 0.67698163 2.010386 747.3633 1.0036847
alpha[4] 1.76067254 3.140392 749.3579 1.0036117
alpha[5] 1.00540867 2.351316 759.6886 1.0040607
alpha[6] 0.84277618 2.190343 744.0685 1.0037104
alpha[7] 0.93667667 2.274858 738.2692 1.0037732
alpha[8] 0.98239879 2.315507 735.8487 1.0034018
alpha[9] 0.06135649 1.409806 771.5404 1.0030400
alpha[10] 0.17629305 1.547322 742.6971 1.0035615
alpha[11] -0.30656376 1.085342 752.0832 1.0036190
beta[1] 0.61306057 1.726937 5533.1425 0.9994772
beta[2] 0.60615024 1.783908 5736.3324 0.9995556
gamma[1] 0.62466093 1.724521 4674.3701 0.9995452
gamma[2] 0.47047709 1.516442 4410.6626 1.0005215
delta[1] 0.66034467 1.647334 5366.8884 0.9994926
delta[2] 0.40148989 1.411885 6169.4977 0.9995869
delta[3] 0.41943426 1.473636 5406.7891 0.9994569
delta[4] 0.59457376 1.613675 5873.7051 0.9997241
rho[1] 0.24758105 1.294939 5933.9098 0.9999407
rho[2] 0.82620151 1.877970 6175.8642 0.9993296
zeta[1] 0.45877262 1.416768 5759.3109 0.9998375
zeta[2] 0.58667190 1.617786 5522.2801 0.9998771
zeta[3] 0.28734899 1.256441 4919.2697 1.0003794
zeta[4] 0.71413244 1.660768 5270.5342 0.9997641Results
Odds ratios
Model A
Phase 1 to Phase 2
ord_mm_tidyA_12 <- ord_mm_fitA_12 %>%
spread_draws(alpha[..], beta[..], c[..]) %>%
mutate(or = exp(beta.1 - beta.2))ggplot(ord_mm_tidyA_12, aes(x = or)) +
geom_histogram() +
geom_vline(xintercept = 1, color = "blue")
table(ord_mm_tidyA_12$or >= 1)/NROW(ord_mm_tidyA_12)
FALSE TRUE
0.542625 0.457375 Phase 1 to Phase 3
ord_mm_tidyA_13 <- ord_mm_fitA_13 %>%
spread_draws(alpha[..], beta[..], c[..]) %>%
mutate(or = exp(beta.1 - beta.2))ggplot(ord_mm_tidyA_13, aes(x = or)) +
geom_histogram() +
geom_vline(xintercept = 1, color = "blue")
table(ord_mm_tidyA_13$or >= 1)/NROW(ord_mm_tidyA_13)
FALSE TRUE
0.00575 0.99425 Model B
Phase 1 to Phase 2
ord_mm_tidyB_12 <- ord_mm_fitB_12 %>%
spread_draws(alpha[..], beta[..], gamma[..], c[..],
delta[..], rho[..], zeta[..]) %>%
mutate(or_men_WA = exp(beta.1 - beta.2 + delta.1 - delta.3 + zeta.1 - zeta.3),
or_women_WA = exp(beta.1 - beta.2 + delta.2 - delta.4 + zeta.1 - zeta.3),
or_men_BHI = exp(beta.1 - beta.2 + delta.1 - delta.3 + zeta.2 - zeta.4),
or_women_BHI = exp(beta.1 - beta.2 + delta.2 - delta.4 + zeta.2 - zeta.4))ggplot(ord_mm_tidyB_12, aes(x = or_men_BHI)) +
geom_histogram() +
geom_vline(xintercept = 1, color = "blue")
table(ord_mm_tidyB_12$or_men_BHI >= 1)/NROW(ord_mm_tidyB_12)
FALSE TRUE
0.7766875 0.2233125 ggplot(ord_mm_tidyB_12, aes(x = or_women_BHI)) +
geom_histogram() +
geom_vline(xintercept = 1, color = "blue")
table(ord_mm_tidyB_12$or_women_BHI >= 1)/NROW(ord_mm_tidyB_12)
FALSE TRUE
0.9945625 0.0054375 ggplot(ord_mm_tidyB_12, aes(x = or_men_WA)) +
geom_histogram() +
geom_vline(xintercept = 1, color = "blue")
table(ord_mm_tidyB_12$or_men_WA >= 1)/NROW(ord_mm_tidyB_12)
FALSE TRUE
0.040875 0.959125 ggplot(ord_mm_tidyB_12, aes(x = or_women_WA)) +
geom_histogram() +
geom_vline(xintercept = 1, color = "blue")
table(ord_mm_tidyB_12$or_women_WA >= 1)/NROW(ord_mm_tidyB_12)
FALSE TRUE
0.8976875 0.1023125 Phase 1 to Phase 3
ord_mm_tidyB_13 <- ord_mm_fitB_13 %>%
spread_draws(alpha[..], beta[..], gamma[..], c[..],
delta[..], rho[..], zeta[..]) %>%
mutate(or_men_WA = exp(beta.1 - beta.2 + delta.1 - delta.3 + zeta.1 - zeta.3),
or_women_WA = exp(beta.1 - beta.2 + delta.2 - delta.4 + zeta.1 - zeta.3),
or_men_BHI = exp(beta.1 - beta.2 + delta.1 - delta.3 + zeta.2 - zeta.4),
or_women_BHI = exp(beta.1 - beta.2 + delta.2 - delta.4 + zeta.2 - zeta.4))ggplot(ord_mm_tidyB_13, aes(x = or_men_BHI)) +
geom_histogram() +
geom_vline(xintercept = 1, color = "blue")
table(ord_mm_tidyB_13$or_men_BHI >= 1)/NROW(ord_mm_tidyB_13)
FALSE TRUE
0.20125 0.79875 ggplot(ord_mm_tidyB_13, aes(x = or_women_BHI)) +
geom_histogram() +
geom_vline(xintercept = 1, color = "blue")
table(ord_mm_tidyB_13$or_women_BHI >= 1)/NROW(ord_mm_tidyB_13)
FALSE TRUE
0.9075 0.0925 ggplot(ord_mm_tidyB_13, aes(x = or_men_WA)) +
geom_histogram() +
geom_vline(xintercept = 1, color = "blue")
table(ord_mm_tidyB_13$or_men_WA >= 1)/NROW(ord_mm_tidyB_13)
FALSE TRUE
0.0005 0.9995 ggplot(ord_mm_tidyB_13, aes(x = or_women_WA)) +
geom_histogram() +
geom_vline(xintercept = 1, color = "blue")
table(ord_mm_tidyB_13$or_women_WA >= 1)/NROW(ord_mm_tidyB_13)
FALSE TRUE
0.573 0.427 Faculty differences for Model A
Phase 1 to Phase 2
ord_mm_tidyA_12 <- ord_mm_fitA_12 %>%
spread_draws(alpha[..], beta[..], c[..]) %>%
mutate(or = exp(beta.1 - beta.2),
A1_1 = 1 - plogis(alpha.1 + beta.1 - c.1),
A2_1 = 1 - plogis(alpha.2 + beta.1 - c.1),
A3_1 = 1 - plogis(alpha.3 + beta.1 - c.1),
A4_1 = 1 - plogis(alpha.4 + beta.1 - c.1),
A5_1 = 1 - plogis(alpha.5 + beta.1 - c.1),
A6_1 = 1 - plogis(alpha.6 + beta.1 - c.1),
A7_1 = 1 - plogis(alpha.7 + beta.1 - c.1),
A8_1 = 1 - plogis(alpha.8 + beta.1 - c.1),
A9_1 = 1 - plogis(alpha.9 + beta.1 - c.1),
A10_1 = 1 - plogis(alpha.10 + beta.1 - c.1),
A11_1 = 1 - plogis(alpha.11 + beta.1 - c.1),
A12_1 = 1 - plogis(alpha.12 + beta.1 - c.1),
A13_1 = 1 - plogis(alpha.13 + beta.1 - c.1),
A14_1 = 1 - plogis(alpha.14 + beta.1 - c.1),
A15_1 = 1 - plogis(alpha.15 + beta.1 - c.1),
A16_1 = 1 - plogis(alpha.16 + beta.1 - c.1),
A17_1 = 1 - plogis(alpha.17 + beta.1 - c.1),
A18_1 = 1 - plogis(alpha.18 + beta.1 - c.1),
A19_1 = 1 - plogis(alpha.19 + beta.1 - c.1),
A20_1 = 1 - plogis(alpha.20 + beta.1 - c.1),
A21_1 = 1 - plogis(alpha.21 + beta.1 - c.1),
A22_1 = 1 - plogis(alpha.22 + beta.1 - c.1),
A23_1 = 1 - plogis(alpha.23 + beta.1 - c.1),
A24_1 = 1 - plogis(alpha.24 + beta.1 - c.1),
A25_1 = 1 - plogis(alpha.25 + beta.1 - c.1),
A1_2 = 1 - plogis(alpha.1 + beta.2 - c.1),
A2_2 = 1 - plogis(alpha.2 + beta.2 - c.1),
A3_2 = 1 - plogis(alpha.3 + beta.2 - c.1),
A4_2 = 1 - plogis(alpha.4 + beta.2 - c.1),
A5_2 = 1 - plogis(alpha.5 + beta.2 - c.1),
A6_2 = 1 - plogis(alpha.6 + beta.2 - c.1),
A7_2 = 1 - plogis(alpha.7 + beta.2 - c.1),
A8_2 = 1 - plogis(alpha.8 + beta.2 - c.1),
A9_2 = 1 - plogis(alpha.9 + beta.2 - c.1),
A10_2 = 1 - plogis(alpha.10 + beta.2 - c.1),
A11_2 = 1 - plogis(alpha.11 + beta.2 - c.1),
A12_2 = 1 - plogis(alpha.12 + beta.2 - c.1),
A13_2 = 1 - plogis(alpha.13 + beta.2 - c.1),
A14_2 = 1 - plogis(alpha.14 + beta.2 - c.1),
A15_2 = 1 - plogis(alpha.15 + beta.2 - c.1),
A16_2 = 1 - plogis(alpha.16 + beta.2 - c.1),
A17_2 = 1 - plogis(alpha.17 + beta.2 - c.1),
A18_2 = 1 - plogis(alpha.18 + beta.2 - c.1),
A19_2 = 1 - plogis(alpha.19 + beta.2 - c.1),
A20_2 = 1 - plogis(alpha.20 + beta.2 - c.1),
A21_2 = 1 - plogis(alpha.21 + beta.2 - c.1),
A22_2 = 1 - plogis(alpha.22 + beta.2 - c.1),
A23_2 = 1 - plogis(alpha.23 + beta.2 - c.1),
A24_2 = 1 - plogis(alpha.24 + beta.2 - c.1),
A25_2 = 1 - plogis(alpha.25 + beta.2 - c.1),
B1_1 = plogis(alpha.1 + beta.1 - c.1) - plogis(alpha.1 + beta.1 - c.2),
B2_1 = plogis(alpha.2 + beta.1 - c.1) - plogis(alpha.2 + beta.1 - c.2),
B3_1 = plogis(alpha.3 + beta.1 - c.1) - plogis(alpha.3 + beta.1 - c.2),
B4_1 = plogis(alpha.4 + beta.1 - c.1) - plogis(alpha.4 + beta.1 - c.2),
B5_1 = plogis(alpha.5 + beta.1 - c.1) - plogis(alpha.5 + beta.1 - c.2),
B6_1 = plogis(alpha.6 + beta.1 - c.1) - plogis(alpha.6 + beta.1 - c.2),
B7_1 = plogis(alpha.7 + beta.1 - c.1) - plogis(alpha.7 + beta.1 - c.2),
B8_1 = plogis(alpha.8 + beta.1 - c.1) - plogis(alpha.8 + beta.1 - c.2),
B9_1 = plogis(alpha.9 + beta.1 - c.1) - plogis(alpha.9 + beta.1 - c.2),
B10_1 = plogis(alpha.10 + beta.1 - c.1) - plogis(alpha.10 + beta.1 - c.2),
B11_1 = plogis(alpha.11 + beta.1 - c.1) - plogis(alpha.11 + beta.1 - c.2),
B12_1 = plogis(alpha.12 + beta.1 - c.1) - plogis(alpha.12 + beta.1 - c.2),
B13_1 = plogis(alpha.13 + beta.1 - c.1) - plogis(alpha.13 + beta.1 - c.2),
B14_1 = plogis(alpha.14 + beta.1 - c.1) - plogis(alpha.14 + beta.1 - c.2),
B15_1 = plogis(alpha.15 + beta.1 - c.1) - plogis(alpha.15 + beta.1 - c.2),
B16_1 = plogis(alpha.16 + beta.1 - c.1) - plogis(alpha.16 + beta.1 - c.2),
B17_1 = plogis(alpha.17 + beta.1 - c.1) - plogis(alpha.17 + beta.1 - c.2),
B18_1 = plogis(alpha.18 + beta.1 - c.1) - plogis(alpha.18 + beta.1 - c.2),
B19_1 = plogis(alpha.19 + beta.1 - c.1) - plogis(alpha.19 + beta.1 - c.2),
B20_1 = plogis(alpha.20 + beta.1 - c.1) - plogis(alpha.20 + beta.1 - c.2),
B21_1 = plogis(alpha.21 + beta.1 - c.1) - plogis(alpha.21 + beta.1 - c.2),
B22_1 = plogis(alpha.22 + beta.1 - c.1) - plogis(alpha.22 + beta.1 - c.2),
B23_1 = plogis(alpha.23 + beta.1 - c.1) - plogis(alpha.23 + beta.1 - c.2),
B24_1 = plogis(alpha.24 + beta.1 - c.1) - plogis(alpha.24 + beta.1 - c.2),
B25_1 = plogis(alpha.25 + beta.1 - c.1) - plogis(alpha.25 + beta.1 - c.2),
B1_2 = plogis(alpha.1 + beta.2 - c.1) - plogis(alpha.1 + beta.2 - c.2),
B2_2 = plogis(alpha.2 + beta.2 - c.1) - plogis(alpha.2 + beta.2 - c.2),
B3_2 = plogis(alpha.3 + beta.2 - c.1) - plogis(alpha.3 + beta.2 - c.2),
B4_2 = plogis(alpha.4 + beta.2 - c.1) - plogis(alpha.4 + beta.2 - c.2),
B5_2 = plogis(alpha.5 + beta.2 - c.1) - plogis(alpha.5 + beta.2 - c.2),
B6_2 = plogis(alpha.6 + beta.2 - c.1) - plogis(alpha.6 + beta.2 - c.2),
B7_2 = plogis(alpha.7 + beta.2 - c.1) - plogis(alpha.7 + beta.2 - c.2),
B8_2 = plogis(alpha.8 + beta.2 - c.1) - plogis(alpha.8 + beta.2 - c.2),
B9_2 = plogis(alpha.9 + beta.2 - c.1) - plogis(alpha.9 + beta.2 - c.2),
B10_2 = plogis(alpha.10 + beta.2 - c.1) - plogis(alpha.10 + beta.2 - c.2),
B11_2 = plogis(alpha.11 + beta.2 - c.1) - plogis(alpha.11 + beta.2 - c.2),
B12_2 = plogis(alpha.12 + beta.2 - c.1) - plogis(alpha.12 + beta.2 - c.2),
B13_2 = plogis(alpha.13 + beta.2 - c.1) - plogis(alpha.13 + beta.2 - c.2),
B14_2 = plogis(alpha.14 + beta.2 - c.1) - plogis(alpha.14 + beta.2 - c.2),
B15_2 = plogis(alpha.15 + beta.2 - c.1) - plogis(alpha.15 + beta.2 - c.2),
B16_2 = plogis(alpha.16 + beta.2 - c.1) - plogis(alpha.16 + beta.2 - c.2),
B17_2 = plogis(alpha.17 + beta.2 - c.1) - plogis(alpha.17 + beta.2 - c.2),
B18_2 = plogis(alpha.18 + beta.2 - c.1) - plogis(alpha.18 + beta.2 - c.2),
B19_2 = plogis(alpha.19 + beta.2 - c.1) - plogis(alpha.19 + beta.2 - c.2),
B20_2 = plogis(alpha.20 + beta.2 - c.1) - plogis(alpha.20 + beta.2 - c.2),
B21_2 = plogis(alpha.21 + beta.2 - c.1) - plogis(alpha.21 + beta.2 - c.2),
B22_2 = plogis(alpha.22 + beta.2 - c.1) - plogis(alpha.22 + beta.2 - c.2),
B23_2 = plogis(alpha.23 + beta.2 - c.1) - plogis(alpha.23 + beta.2 - c.2),
B24_2 = plogis(alpha.24 + beta.2 - c.1) - plogis(alpha.24 + beta.2 - c.2),
B25_2 = plogis(alpha.25 + beta.2 - c.1) - plogis(alpha.25 + beta.2 - c.2),
C1_1 = plogis(alpha.1 + beta.1 - c.2) - plogis(alpha.1 + beta.1 - c.3),
C2_1 = plogis(alpha.2 + beta.1 - c.2) - plogis(alpha.2 + beta.1 - c.3),
C3_1 = plogis(alpha.3 + beta.1 - c.2) - plogis(alpha.3 + beta.1 - c.3),
C4_1 = plogis(alpha.4 + beta.1 - c.2) - plogis(alpha.4 + beta.1 - c.3),
C5_1 = plogis(alpha.5 + beta.1 - c.2) - plogis(alpha.5 + beta.1 - c.3),
C6_1 = plogis(alpha.6 + beta.1 - c.2) - plogis(alpha.6 + beta.1 - c.3),
C7_1 = plogis(alpha.7 + beta.1 - c.2) - plogis(alpha.7 + beta.1 - c.3),
C8_1 = plogis(alpha.8 + beta.1 - c.2) - plogis(alpha.8 + beta.1 - c.3),
C9_1 = plogis(alpha.9 + beta.1 - c.2) - plogis(alpha.9 + beta.1 - c.3),
C10_1 = plogis(alpha.10 + beta.1 - c.2) - plogis(alpha.10 + beta.1 - c.3),
C11_1 = plogis(alpha.11 + beta.1 - c.2) - plogis(alpha.11 + beta.1 - c.3),
C12_1 = plogis(alpha.12 + beta.1 - c.2) - plogis(alpha.12 + beta.1 - c.3),
C13_1 = plogis(alpha.13 + beta.1 - c.2) - plogis(alpha.13 + beta.1 - c.3),
C14_1 = plogis(alpha.14 + beta.1 - c.2) - plogis(alpha.14 + beta.1 - c.3),
C15_1 = plogis(alpha.15 + beta.1 - c.2) - plogis(alpha.15 + beta.1 - c.3),
C16_1 = plogis(alpha.16 + beta.1 - c.2) - plogis(alpha.16 + beta.1 - c.3),
C17_1 = plogis(alpha.17 + beta.1 - c.2) - plogis(alpha.17 + beta.1 - c.3),
C18_1 = plogis(alpha.18 + beta.1 - c.2) - plogis(alpha.18 + beta.1 - c.3),
C19_1 = plogis(alpha.19 + beta.1 - c.2) - plogis(alpha.19 + beta.1 - c.3),
C20_1 = plogis(alpha.20 + beta.1 - c.2) - plogis(alpha.20 + beta.1 - c.3),
C21_1 = plogis(alpha.21 + beta.1 - c.2) - plogis(alpha.21 + beta.1 - c.3),
C22_1 = plogis(alpha.22 + beta.1 - c.2) - plogis(alpha.22 + beta.1 - c.3),
C23_1 = plogis(alpha.23 + beta.1 - c.2) - plogis(alpha.23 + beta.1 - c.3),
C24_1 = plogis(alpha.24 + beta.1 - c.2) - plogis(alpha.24 + beta.1 - c.3),
C25_1 = plogis(alpha.25 + beta.1 - c.2) - plogis(alpha.25 + beta.1 - c.3),
C1_2 = plogis(alpha.1 + beta.2 - c.2) - plogis(alpha.1 + beta.2 - c.3),
C2_2 = plogis(alpha.2 + beta.2 - c.2) - plogis(alpha.2 + beta.2 - c.3),
C3_2 = plogis(alpha.3 + beta.2 - c.2) - plogis(alpha.3 + beta.2 - c.3),
C4_2 = plogis(alpha.4 + beta.2 - c.2) - plogis(alpha.4 + beta.2 - c.3),
C5_2 = plogis(alpha.5 + beta.2 - c.2) - plogis(alpha.5 + beta.2 - c.3),
C6_2 = plogis(alpha.6 + beta.2 - c.2) - plogis(alpha.6 + beta.2 - c.3),
C7_2 = plogis(alpha.7 + beta.2 - c.2) - plogis(alpha.7 + beta.2 - c.3),
C8_2 = plogis(alpha.8 + beta.2 - c.2) - plogis(alpha.8 + beta.2 - c.3),
C9_2 = plogis(alpha.9 + beta.2 - c.2) - plogis(alpha.9 + beta.2 - c.3),
C10_2 = plogis(alpha.10 + beta.2 - c.2) - plogis(alpha.10 + beta.2 - c.3),
C11_2 = plogis(alpha.11 + beta.2 - c.2) - plogis(alpha.11 + beta.2 - c.3),
C12_2 = plogis(alpha.12 + beta.2 - c.2) - plogis(alpha.12 + beta.2 - c.3),
C13_2 = plogis(alpha.13 + beta.2 - c.2) - plogis(alpha.13 + beta.2 - c.3),
C14_2 = plogis(alpha.14 + beta.2 - c.2) - plogis(alpha.14 + beta.2 - c.3),
C15_2 = plogis(alpha.15 + beta.2 - c.2) - plogis(alpha.15 + beta.2 - c.3),
C16_2 = plogis(alpha.16 + beta.2 - c.2) - plogis(alpha.16 + beta.2 - c.3),
C17_2 = plogis(alpha.17 + beta.2 - c.2) - plogis(alpha.17 + beta.2 - c.3),
C18_2 = plogis(alpha.18 + beta.2 - c.2) - plogis(alpha.18 + beta.2 - c.3),
C19_2 = plogis(alpha.19 + beta.2 - c.2) - plogis(alpha.19 + beta.2 - c.3),
C20_2 = plogis(alpha.20 + beta.2 - c.2) - plogis(alpha.20 + beta.2 - c.3),
C21_2 = plogis(alpha.21 + beta.2 - c.2) - plogis(alpha.21 + beta.2 - c.3),
C22_2 = plogis(alpha.22 + beta.2 - c.2) - plogis(alpha.22 + beta.2 - c.3),
C23_2 = plogis(alpha.23 + beta.2 - c.2) - plogis(alpha.23 + beta.2 - c.3),
C24_2 = plogis(alpha.24 + beta.2 - c.2) - plogis(alpha.24 + beta.2 - c.3),
C25_2 = plogis(alpha.25 + beta.2 - c.2) - plogis(alpha.25 + beta.2 - c.3),
DFW1_1 = plogis(alpha.1 + beta.1 - c.3),
DFW2_1 = plogis(alpha.2 + beta.1 - c.3),
DFW3_1 = plogis(alpha.3 + beta.1 - c.3),
DFW4_1 = plogis(alpha.4 + beta.1 - c.3),
DFW5_1 = plogis(alpha.5 + beta.1 - c.3),
DFW6_1 = plogis(alpha.6 + beta.1 - c.3),
DFW7_1 = plogis(alpha.7 + beta.1 - c.3),
DFW8_1 = plogis(alpha.8 + beta.1 - c.3),
DFW9_1 = plogis(alpha.9 + beta.1 - c.3),
DFW10_1 = plogis(alpha.10 + beta.1 - c.3),
DFW11_1 = plogis(alpha.11 + beta.1 - c.3),
DFW12_1 = plogis(alpha.12 + beta.1 - c.3),
DFW13_1 = plogis(alpha.13 + beta.1 - c.3),
DFW14_1 = plogis(alpha.14 + beta.1 - c.3),
DFW15_1 = plogis(alpha.15 + beta.1 - c.3),
DFW16_1 = plogis(alpha.16 + beta.1 - c.3),
DFW17_1 = plogis(alpha.17 + beta.1 - c.3),
DFW18_1 = plogis(alpha.18 + beta.1 - c.3),
DFW19_1 = plogis(alpha.19 + beta.1 - c.3),
DFW20_1 = plogis(alpha.20 + beta.1 - c.3),
DFW21_1 = plogis(alpha.21 + beta.1 - c.3),
DFW22_1 = plogis(alpha.22 + beta.1 - c.3),
DFW23_1 = plogis(alpha.23 + beta.1 - c.3),
DFW24_1 = plogis(alpha.24 + beta.1 - c.3),
DFW25_1 = plogis(alpha.25 + beta.1 - c.3),
DFW1_2 = plogis(alpha.1 + beta.2 - c.3),
DFW2_2 = plogis(alpha.2 + beta.2 - c.3),
DFW3_2 = plogis(alpha.3 + beta.2 - c.3),
DFW4_2 = plogis(alpha.4 + beta.2 - c.3),
DFW5_2 = plogis(alpha.5 + beta.2 - c.3),
DFW6_2 = plogis(alpha.6 + beta.2 - c.3),
DFW7_2 = plogis(alpha.7 + beta.2 - c.3),
DFW8_2 = plogis(alpha.8 + beta.2 - c.3),
DFW9_2 = plogis(alpha.9 + beta.2 - c.3),
DFW10_2 = plogis(alpha.10 + beta.2 - c.3),
DFW11_2 = plogis(alpha.11 + beta.2 - c.3),
DFW12_2 = plogis(alpha.12 + beta.2 - c.3),
DFW13_2 = plogis(alpha.13 + beta.2 - c.3),
DFW14_2 = plogis(alpha.14 + beta.2 - c.3),
DFW15_2 = plogis(alpha.15 + beta.2 - c.3),
DFW16_2 = plogis(alpha.16 + beta.2 - c.3),
DFW17_2 = plogis(alpha.17 + beta.2 - c.3),
DFW18_2 = plogis(alpha.18 + beta.2 - c.3),
DFW19_2 = plogis(alpha.19 + beta.2 - c.3),
DFW20_2 = plogis(alpha.20 + beta.2 - c.3),
DFW21_2 = plogis(alpha.21 + beta.2 - c.3),
DFW22_2 = plogis(alpha.22 + beta.2 - c.3),
DFW23_2 = plogis(alpha.23 + beta.2 - c.3),
DFW24_2 = plogis(alpha.24 + beta.2 - c.3),
DFW25_2 = plogis(alpha.25 + beta.2 - c.3))
ord_mm_tidyA_12ord_mm_tidy_diffA_12 <- ord_mm_tidyA_12 %>%
transmute(A1_diff = A1_2 - A1_1,
A2_diff = A2_2 - A2_1,
A3_diff = A3_2 - A3_1,
A4_diff = A4_2 - A4_1,
A5_diff = A5_2 - A5_1,
A6_diff = A6_2 - A6_1,
A7_diff = A7_2 - A7_1,
A8_diff = A8_2 - A8_1,
A9_diff = A9_2 - A9_1,
A10_diff = A10_2 - A10_1,
A11_diff = A11_2 - A11_1,
A12_diff = A12_2 - A12_1,
A13_diff = A13_2 - A13_1,
A14_diff = A14_2 - A14_1,
A15_diff = A15_2 - A15_1,
A16_diff = A16_2 - A16_1,
A17_diff = A17_2 - A17_1,
A18_diff = A18_2 - A18_1,
A19_diff = A19_2 - A19_1,
A20_diff = A20_2 - A20_1,
A21_diff = A21_2 - A21_1,
A22_diff = A22_2 - A22_1,
A23_diff = A23_2 - A23_1,
A24_diff = A24_2 - A24_1,
A25_diff = A25_2 - A25_1,
B1_diff = B1_2 - B1_1,
B2_diff = B2_2 - B2_1,
B3_diff = B3_2 - B3_1,
B4_diff = B4_2 - B4_1,
B5_diff = B5_2 - B5_1,
B6_diff = B6_2 - B6_1,
B7_diff = B7_2 - B7_1,
B8_diff = B8_2 - B8_1,
B9_diff = B9_2 - B9_1,
B10_diff = B10_2 - B10_1,
B11_diff = B11_2 - B11_1,
B12_diff = B12_2 - B12_1,
B13_diff = B13_2 - B13_1,
B14_diff = B14_2 - B14_1,
B15_diff = B15_2 - B15_1,
B16_diff = B16_2 - B16_1,
B17_diff = B17_2 - B17_1,
B18_diff = B18_2 - B18_1,
B19_diff = B19_2 - B19_1,
B20_diff = B20_2 - B20_1,
B21_diff = B21_2 - B21_1,
B22_diff = B22_2 - B22_1,
B23_diff = B23_2 - B23_1,
B24_diff = B24_2 - B24_1,
B25_diff = B25_2 - B25_1,
C1_diff = C1_2 - C1_1,
C2_diff = C2_2 - C2_1,
C3_diff = C3_2 - C3_1,
C4_diff = C4_2 - C4_1,
C5_diff = C5_2 - C5_1,
C6_diff = C6_2 - C6_1,
C7_diff = C7_2 - C7_1,
C8_diff = C8_2 - C8_1,
C9_diff = C9_2 - C9_1,
C10_diff = C10_2 - C10_1,
C11_diff = C11_2 - C11_1,
C12_diff = C12_2 - C12_1,
C13_diff = C13_2 - C13_1,
C14_diff = C14_2 - C14_1,
C15_diff = C15_2 - C15_1,
C16_diff = C16_2 - C16_1,
C17_diff = C17_2 - C17_1,
C18_diff = C18_2 - C18_1,
C19_diff = C19_2 - C19_1,
C20_diff = C20_2 - C20_1,
C21_diff = C21_2 - C21_1,
C22_diff = C22_2 - C22_1,
C23_diff = C23_2 - C23_1,
C24_diff = C24_2 - C24_1,
C25_diff = C25_2 - C25_1,
DFW1_diff = DFW1_2 - DFW1_1,
DFW2_diff = DFW2_2 - DFW2_1,
DFW3_diff = DFW3_2 - DFW3_1,
DFW4_diff = DFW4_2 - DFW4_1,
DFW5_diff = DFW5_2 - DFW5_1,
DFW6_diff = DFW6_2 - DFW6_1,
DFW7_diff = DFW7_2 - DFW7_1,
DFW8_diff = DFW8_2 - DFW8_1,
DFW9_diff = DFW9_2 - DFW9_1,
DFW10_diff = DFW10_2 - DFW10_1,
DFW11_diff = DFW11_2 - DFW11_1,
DFW12_diff = DFW12_2 - DFW12_1,
DFW13_diff = DFW13_2 - DFW13_1,
DFW14_diff = DFW14_2 - DFW14_1,
DFW15_diff = DFW15_2 - DFW15_1,
DFW16_diff = DFW16_2 - DFW16_1,
DFW17_diff = DFW17_2 - DFW17_1,
DFW18_diff = DFW18_2 - DFW18_1,
DFW19_diff = DFW19_2 - DFW19_1,
DFW20_diff = DFW20_2 - DFW20_1,
DFW21_diff = DFW21_2 - DFW21_1,
DFW22_diff = DFW22_2 - DFW22_1,
DFW23_diff = DFW23_2 - DFW23_1,
DFW24_diff = DFW24_2 - DFW24_1,
DFW25_diff = DFW25_2 - DFW25_1)
ord_mm_tidy_diffA_12